#!pip install mca
#!pip install prince
import pandas as pd
import numpy as np
#from sklearn.cluster import KMeans, AgglomerativeClustering
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.preprocessing import MinMaxScaler
from sklearn.decomposition import PCA
import mca
import prince
import itertools
import scipy.stats as ss
Jak dowiadujemy się z opisu danych są one obrobioną nieco wersją zbioru USCensus1990raw. Zawierał on 1 % próbki z pełnego spisu powszechnego z 1990r. W naszym pliku rekordy zostały losowo przemieszane. Aby uzyskać zbiór, którym się zajmujemy zostały wykonane następujące operacje:
Można więc powiedzieć, że na danych został już wykonany preprocessing.
Pierwszym atrybutem jest caseid i powinnien on być ignorowany podczas analizy.
Wszystkie dane są numeryczne.
X = pd.read_csv('USCensus1990.data.txt')
X.head()
| caseid | dAge | dAncstry1 | dAncstry2 | iAvail | iCitizen | iClass | dDepart | iDisabl1 | iDisabl2 | ... | iTmpabsnt | dTravtime | iVietnam | dWeek89 | iWork89 | iWorklwk | iWWII | iYearsch | iYearwrk | dYrsserv | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 10000 | 5 | 0 | 1 | 0 | 0 | 5 | 3 | 2 | 2 | ... | 0 | 5 | 0 | 2 | 1 | 1 | 0 | 11 | 1 | 0 |
| 1 | 10001 | 6 | 1 | 1 | 0 | 0 | 7 | 5 | 2 | 2 | ... | 0 | 1 | 0 | 2 | 1 | 1 | 0 | 5 | 1 | 0 |
| 2 | 10002 | 3 | 1 | 2 | 0 | 0 | 7 | 4 | 2 | 2 | ... | 0 | 2 | 0 | 2 | 1 | 1 | 0 | 10 | 1 | 0 |
| 3 | 10003 | 4 | 1 | 2 | 0 | 0 | 1 | 3 | 2 | 2 | ... | 0 | 1 | 0 | 1 | 1 | 1 | 0 | 10 | 1 | 0 |
| 4 | 10004 | 7 | 1 | 1 | 0 | 0 | 0 | 0 | 2 | 2 | ... | 3 | 0 | 0 | 0 | 2 | 2 | 0 | 5 | 6 | 0 |
5 rows × 69 columns
X.shape
(2458285, 69)
Sprawdźmy czy w zbiorze są jakieś braki danych:
X.isna().sum().sum()
0
X.isnull().sum().sum()
0
Wygląda na to, że w naszym zbiorze nie ma braków danych.
X.duplicated().sum()
0
Nie ma też zduplikowanych wierszy.
Sprawdzimy jeszcze jakie są unikalne wartości w każdej z kolumn:
pd.set_option('display.max_rows', None)
print(X.apply(lambda col: col.unique()))
caseid [10000, 10001, 10002, 10003, 10004, 10005, 100... dAge [5, 6, 3, 4, 7, 1, 2, 0] dAncstry1 [0, 1, 11, 3, 2, 9, 8, 10, 4, 6, 7, 5] dAncstry2 [1, 2, 12, 3, 4, 10, 7, 5, 11, 8, 6, 9] iAvail [0, 3, 4, 2, 1] iCitizen [0, 3, 4, 2, 1] iClass [5, 7, 1, 0, 6, 3, 4, 2, 8, 9] dDepart [3, 5, 4, 0, 2, 1] iDisabl1 [2, 0, 1] iDisabl2 [2, 0, 1] iEnglish [1, 0, 2, 3, 4] iFeb55 [0, 1] iFertil [1, 3, 0, 4, 7, 2, 5, 6, 12, 8, 13, 10, 9, 11] dHispanic [0, 1, 9, 2, 6, 3, 4, 8, 7, 5] dHour89 [4, 1, 3, 0, 5, 2] dHours [3, 1, 4, 0, 5, 2] iImmigr [0, 5, 1, 8, 3, 10, 6, 7, 4, 2, 9] dIncome1 [2, 1, 0, 3, 4] dIncome2 [0, 1] dIncome3 [0, 1] dIncome4 [1, 0] dIncome5 [0, 1] dIncome6 [0, 1] dIncome7 [0, 1] dIncome8 [0, 1] dIndustry [10, 4, 1, 0, 9, 7, 6, 3, 5, 8, 11, 12, 2] iKorean [0, 1] iLang1 [1, 2, 0] iLooking [0, 2, 1] iMarital [1, 0, 4, 2, 3] iMay75880 [0, 1] iMeans [1, 10, 0, 11, 2, 4, 8, 5, 12, 3, 7, 9, 6] iMilitary [4, 0, 2, 3, 1] iMobility [2, 1, 0] iMobillim [2, 0, 1] dOccup [3, 2, 4, 0, 6, 1, 5, 8, 7] iOthrserv [0, 1] iPerscare [2, 0, 1] dPOB [0, 3, 4, 6, 2, 1, 5] dPoverty [2, 1, 0] dPwgt1 [1, 2, 0, 3] iRagechld [4, 2, 0, 3, 1] dRearning [3, 2, 0, 4, 5, 1] iRelat1 [0, 1, 2, 10, 12, 9, 7, 4, 3, 13, 11, 6, 8, 5] iRelat2 [0, 6, 9, 1, 8, 2, 4, 3, 7, 5] iRemplpar [0, 121, 112, 211, 111, 122, 222, 221, 223, 11... iRiders [3, 1, 0, 2, 4, 8, 6, 5, 7] iRlabor [1, 6, 0, 3, 4, 2, 5] iRownchld [0, 1] dRpincome [3, 2, 0, 4, 5, 1] iRPOB [22, 10, 21, 23, 51, 24, 52, 40, 31, 35, 36, 3... iRrelchld [0, 1] iRspouse [3, 1, 6, 4, 0, 5, 2] iRvetserv [0, 1, 3, 8, 10, 4, 7, 2, 5, 9, 6, 11] iSchool [1, 2, 3, 0] iSept80 [0, 1] iSex [1, 0] iSubfam1 [0, 2, 3, 1] iSubfam2 [0, 1, 2, 3] iTmpabsnt [0, 3, 1, 2] dTravtime [5, 1, 2, 0, 4, 3, 6] iVietnam [0, 1] dWeek89 [2, 1, 0] iWork89 [1, 2, 0] iWorklwk [1, 2, 0] iWWII [0, 1] iYearsch [11, 5, 10, 4, 8, 1, 7, 6, 12, 14, 15, 2, 0, 9... iYearwrk [1, 6, 0, 7, 4, 2, 3, 5] dYrsserv [0, 1, 2] dtype: object
pd.set_option('display.max_rows', 10) #back to default
Oryginalny zbiór jest duży, ograniczymy go do losowej jednoprocentowej próbki, dodatkowo usuniemy zmienną caseid.
df = X.sample(frac=0.01, random_state=17)
df.drop('caseid', inplace = True, axis = 1)
df.shape
(24583, 68)
Przejdziemy teraz do głębszej analizy zmiennych zawartych w naszej ramce. Przyjrzymy sie mapie korelacji a także zależnościom między zmiennymi.
plt.figure(figsize=(40, 40))
sns.heatmap(df.corr(method = 'pearson'), annot=True)
plt.show()
Z powyższej mapy możemy odczytać, że w ramce danych występują wysoko skorelowane zmienne. Podzbędziemy się tych o najwyższym module wartości korelacji (czyli większym od 0.7).
correlated_features = set()
corr = df.corr(method = 'spearman')
for i in range(len(corr .columns)):
for j in range(i):
if abs(corr.iloc[i, j]) >= 0.70:
colname = corr.columns[i]
correlated_features.add(colname)
print(len(correlated_features))
print(correlated_features)
29
{'iMeans', 'iYearwrk', 'dPOB', 'iSubfam2', 'iSex', 'iRownchld', 'iRiders', 'iDisabl2', 'dRearning', 'iYearsch', 'iTmpabsnt', 'iRemplpar', 'dIncome1', 'dOccup', 'dTravtime', 'dHours', 'iRlabor', 'dYrsserv', 'dRpincome', 'iMilitary', 'dWeek89', 'iWork89', 'dHour89', 'iPerscare', 'iWorklwk', 'iImmigr', 'iRrelchld', 'dIndustry', 'iMobillim'}
df.drop(labels=correlated_features, axis=1, inplace=True)
df.shape
(24583, 39)
Po obróbce zostało 38 kolumn w ramce. Jako, że wszystkie z nich są kategoryczne, zbadamy teraz ich rozkłady na grupowym histogramie.
df1 = X.sample(frac=0.01, random_state=17)
df1.drop('caseid', inplace = True, axis = 1)
df_mca = mca.MCA(df1)
Sprawdźmy najpierw do jakiej liczby zmiennych możemy zredukować:
plt.plot(np.cumsum(df_mca.expl_var(greenacre=False, N=68)))
plt.xlabel('number of components')
plt.ylabel('cumulative explained variance');
Z wykresu wynika, że najbardziej odpowiednią wartością będzie 10.
mca = prince.MCA(n_components=10)
mca = mca.fit(df1)
mca = mca.transform(df1)
mca
| 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | |
|---|---|---|---|---|---|---|---|---|---|---|
| 1551174 | -0.695760 | -0.366005 | -0.118154 | -0.232889 | -0.092356 | 0.174044 | 0.000954 | 0.087259 | -0.093411 | -0.119080 |
| 1487409 | 1.077608 | -0.013721 | -0.039326 | -0.032668 | -0.014635 | 0.063703 | 0.053992 | -0.033713 | -0.036856 | -0.030213 |
| 436156 | -0.648278 | 0.923192 | -0.427529 | -0.059384 | -0.139972 | 0.159919 | -0.012562 | 0.067955 | -0.045496 | -0.023146 |
| 1088689 | -0.695765 | -0.460353 | 0.121820 | -0.019130 | -0.067644 | 0.107694 | 0.007207 | 0.031381 | -0.134219 | -0.017786 |
| 678665 | -0.561437 | 0.819894 | -0.151287 | -0.156990 | -0.084995 | -0.053536 | -0.403601 | 0.055249 | 0.300525 | 0.011119 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 77165 | -0.651350 | 0.868671 | -0.389745 | -0.298434 | -0.183215 | 0.337137 | 0.045245 | 0.046961 | -0.005184 | -0.022022 |
| 823271 | -0.615616 | 0.668123 | -0.125566 | -0.235723 | -0.191776 | 0.367406 | 0.012391 | -0.011987 | 0.047139 | -0.065916 |
| 1797933 | -0.651221 | 0.496347 | -0.106328 | -0.112082 | -0.045619 | -0.654337 | 0.081591 | 0.146347 | -0.834846 | 0.908348 |
| 1996241 | -0.679699 | -0.315659 | 0.067654 | -0.178622 | -0.117867 | 0.270970 | -0.043009 | 0.055422 | -0.072660 | -0.099030 |
| 479824 | -0.642057 | -0.282944 | -0.163976 | 0.558653 | 0.347452 | -0.335598 | 0.242525 | -0.247698 | -0.033250 | -0.432841 |
24583 rows × 10 columns
def cramers_v(x, y):
confusion_matrix = pd.crosstab(x,y)
chi2 = ss.chi2_contingency(confusion_matrix)[0]
n = confusion_matrix.sum().sum()
phi2 = chi2/n
r,k = confusion_matrix.shape
phi2corr = max(0, phi2-((k-1)*(r-1))/(n-1))
rcorr = r-((r-1)**2)/(n-1)
kcorr = k-((k-1)**2)/(n-1)
return np.sqrt(phi2corr/min((kcorr-1),(rcorr-1)))
correlated_features = set()
for i in range(len(df1.columns)):
for j in range(i):
if cramers_v(df1[df1.columns[i]], df1[df1.columns[j]]) >= 0.70:
colname = df1.columns[i]
correlated_features.add(colname)
print(correlated_features)
{'iYearwrk', 'iRvetserv', 'dPOB', 'iSex', 'iWWII', 'iMobility', 'iRownchld', 'iVietnam', 'iDisabl2', 'dRearning', 'iYearsch', 'iTmpabsnt', 'iRspouse', 'dIncome5', 'dIncome2', 'iSept80', 'iRlabor', 'dYrsserv', 'dRpincome', 'iMilitary', 'dWeek89', 'iWork89', 'iPerscare', 'iLang1', 'iWorklwk', 'iRPOB', 'iSchool', 'iRrelchld', 'iMobillim'}
cols = df1.columns.values.tolist()
corrM = np.zeros((len(cols),len(cols)))
# there's probably a nice pandas way to do this
for col1, col2 in itertools.combinations(cols, 2):
idx1, idx2 = cols.index(col1), cols.index(col2)
corrM[idx1, idx2] = cramers_v(df1[col1], df1[col2])
corrM[idx2, idx1] = corrM[idx1, idx2]
corr = pd.DataFrame(corrM, index=cols, columns=cols)
fig, ax = plt.subplots(figsize=(40, 40))
ax = sns.heatmap(corr, annot=True, ax=ax); ax.set_title("Cramer V Correlation between Variables");
df.hist(figsize=(15,20))
plt.show()
Możemy zauważyć, że większość zmiennych ma bardzo nierównomierny rozkład, w którym przeważa jedna lub dwie wartości. Jednynie zmienna dAge przypomina rozkład jednostajny.
Następnie przyjrzemy się wykresom wizualizującym relacje między zmiennymi, aby odkryć jakieś zależności, które mogą nam pomóc przy doborze klastrów. Przyjrzemy się bubble plotom, które pokazują liczności przecięć dwóch zmiennych.
def HeatMap(column1, column2):
df1 = df.groupby([column1, column2]).size().reset_index(name='counts')
df1 = pd.pivot(df1, index = 'race', columns = 'sex', values = 'counts')
sns.heatmap(df1, annot=True, cmap='Reds', ax=ax1, fmt='g', norm=LogNorm())
def BubblePlot(column1, column2):
df1 = df.groupby([column1, column2]).size().reset_index(name='counts')
plt.figure(figsize = (12, 8))
plt.scatter(df1[column1], df1[column2], s=df1['counts'], alpha = 0.7)
plt.xlabel(column1)
plt.gca().xaxis.set_major_locator(plt.MultipleLocator(1))
plt.gca().yaxis.set_major_locator(plt.MultipleLocator(1))
plt.ylabel(column2)
df1 = df.groupby(['iMarital', 'iClass']).size().reset_index(name='counts')
BubblePlot('iMarital', 'iClass')
Legenda:
| iMarital (marital status) | iClass (class of worker) | |
|---|---|---|
| 0 | Now Married, Except Separated | N/a Less Than 16 Yrs. Old/unemp. Who Never Worked |
| 1 | Widowed | Emp. of a Private for Profit Company |
| 2 | Divorced | Emp. of a Private Not for Profit |
| 3 | Separated | Local Gov. Emp. City, County, Etc. |
| 4 | Never Married or Under 15 Yrs. Old | State Gov. Emp. |
| 5 | Federal Gov. Emp. | |
| 6 | Self Emp. in Own Not Incorp.d Business | |
| 7 | Self Emp. in Own Incorp.d Business | |
| 8 | Working Without Pay in Family Business or Farm | |
| 9 | Unemp., Last Worked in 1984 or Earlier |
Widać, że przeważają grupy osób, które są w związku małżeńskim (1 na osi x) nigdy nie były w sformalizowanym związku lub poniiżej 15 roku życia (4 na osi x). Dodatkowo widać, że w grupie 4 na osi x, jest dużo osób, które nigdy nie były zatrudnione lub mają mniej niż 16 lat (0 na osi y), co jest dosyć intuicyjne. Dodatkowo wśród osób zamężnych jest przewaga osób pracujących prywatnie (1 na osi y)
BubblePlot('dAge', "iRelat1")
Legenda:
| iRelat1 (related or not) | dAge | |
|---|---|---|
| 0 | Householder | 0 |
| 1 | Husband/wife | <13 |
| 2 | Son/daughter | 13-19 |
| 3 | Stepson/stepdaughter | 20-29 |
| 4 | Brother/sister | 30-39 |
| 5 | Father/mother | 40-49 |
| 6 | Grandchild | 50-64 |
| 7 | Other Rel. | >65 |
| 8 | Roomer/boarder/foster Child | |
| 9 | Housemate/roommate | |
| 10 | Unmarried Partner | |
| 11 | Other Nonrel. | |
| 12 | Instit. Person | |
| 13 | Other Pers. in Grp. Qtrs. |
Na wykresie ze zmiennymi dAge i iRelat1 widać niesamowicie ciekawe i oszałamiające rzeczy. Widać, że dorosłe osoby (oznaczone numerami większymi niż 2 na osi x) są głowami rodziny (0 na osi y) i są zakawalifikowani jako mąż lub żona (1 na osi y). Natomiast osoby młodsze (1 na osi x) są kwalifikowani przede wszytskim jako dzieci, co jest spójne.
BubblePlot('iEnglish', 'iSchool')
Legenda:
| iEnglish (ability to speak English) | iSchool (school enrollment) | |
|---|---|---|
| 0 | N/a Less Than 5 Yrs. Old/speaks Only English | N/a Less Than 3 Yrs. Old |
| 1 | Very well | Not Attending School |
| 2 | Well | Yes, Pub. School, Pub. Coll. |
| 3 | Not Well | Yes, Private School, Private Coll. |
| 4 | Not at all |
Największa grupą wśród rozważanych stanową osoby, które mówią tylko po angielsku. Wśród nich najwięcej jest osób, które nie chodzą do szkoły, czyli jak można wnioskować - dorosłych. Można też zauważyć, że więcej osób chodzi do publicznych, niż prywatnych placówek.
Aby znormalizować naszą ramkę danych, postanowiliśmy przeskalować dane, aby zawierały się tylko w zbiorze [0,1]. Zmienne nadal pozostaną dyksretne, ale będzie zdecydowanie łatwiej analizować dobierane klastry w późniejszych krokach. Mamy na uwadze również to, że takie skalowanie zmieni dobór klastrów, bo odległości między punktami będą ewaluowane inaczej.
for col in df.columns:
min_max_scaler1 = MinMaxScaler()
scaled_col1 = min_max_scaler1.fit_transform(df[col].values.reshape(-1, 1))
df[col] = scaled_col1
from sklearn.cluster import KMeans, AgglomerativeClustering
from sklearn.metrics import silhouette_score, davies_bouldin_score
from scipy.cluster import hierarchy
from sklearn import cluster, datasets, mixture
from sklearn_extra.cluster import KMedoids
from yellowbrick.cluster import SilhouetteVisualizer
clus_num = range(1, 20)
scores = []
for k in clus_num:
model = KMeans(n_clusters = k)
model.fit(df)
scores.append(model.score(df) * -1)
plt.plot(clus_num, scores, '-o')
plt.xlabel('number of clusters, k')
plt.ylabel('Within-cluster sum of squares')
plt.xticks(clus_num)
plt.show()
Na początek zbadaliśmy metodą łokcia, jaka liczba klastrów będzie optymalna. Wykres scorów dla liczby klastrów od 1 do 19 jest słabo stromy, więc trudno wybrać jednoznacznie punkt przegięcie. Ustalmy jednak, ze znajduje się on w okolicach wartośći od 5 d 9
Uznaliśmy, ze dobrą metodą na wizualizacje jak model dzieli dane na klastry, jest wykorzystanie centroidów i porównywanie ich w zależności od klastra dla każdej ze zmiennych.
def CentroidsPlot(centroids):
cent_df = pd.DataFrame(columns= ['cluster', 'variable', 'centroid'])
rows = centroids.shape[0]
clusters = [x+1 for x in range(rows)]
for colname in centroids.columns:
data = {'cluster' : clusters,
'variable' : [colname for x in range(rows)],
'centroid' : centroids[colname]}
cent_df = pd.concat([cent_df, pd.DataFrame(data)])
plt.figure(figsize=(18, 8))
scatter = plt.scatter(cent_df.variable, cent_df.centroid, c=cent_df.cluster, cmap='viridis', s = 200, alpha = 0.7)
plt.legend(handles=scatter.legend_elements()[0],
labels=clusters,
title="cluster")
plt.xticks(rotation = 90, fontsize = 15)
Na początek zbudujemy modele KMeans dla znalezionych wcześniej optymalnych liczb klastrów
model = KMeans(n_clusters = 5)
model.fit(df)
labs = model.labels_
centroids = model.cluster_centers_
centroids = pd.DataFrame(columns = df.columns, data = centroids)
CentroidsPlot(centroids)
model = KMeans(n_clusters = 7)
model.fit(df)
labs = model.labels_
centroids = model.cluster_centers_
centroids = pd.DataFrame(columns = df.columns, data = centroids)
CentroidsPlot(centroids)
model = KMeans(n_clusters = 9)
model.fit(df)
labs = model.labels_
centroids = model.cluster_centers_
centroids = pd.DataFrame(columns = df.columns, data = centroids)
CentroidsPlot(centroids)
Widać, że dobór większej liczby klastrów sprawia problemy, ponieważ kolejne klastry są do siebie za bardzo zbliżone. Dla 5 klastrów widać, że grupy mają różne dobrane centroidy, przez co wydaje się, że rzeczywiście dobrze rozdzielają rekory na różne klastry. Na pewno widać, że klaster oznaczony jako 3 dobrze oddziela młodzież i dzieci, natomiast klaster 4 ludzi starszych.